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ABSTRACT 

Bressert et al. recently showed that the surface density distribution of low-mass, young stellar 
objects in the solar neighbourhood is approximately lognormal. The authors conclude that the 
star formation process is hierarchical and that only a small fraction of stars form in dense star 
clusters. Here we show that the peak and the width of the density distribution is also what 
follows if all stars form in bound clusters which are not significantly affected by the presence 
of gas and expand by two-body relaxation. The peak of the surface density distribution is sim- 
ply obtained from the typical ages (few Myrs) and cluster membership number (few hundred) 
typifying nearby star forming regions. This result depends weakly on initial cluster sizes, 
provided that they are sufficiently dense (initial half mass radius of < 0.3 pc) for dynamical 
evolution to be important at an age of a few Myrs. We conclude that the degeneracy of the 
YSO surface density distribution complicates its use as a diagnostic of the stellar formation 
environment. 
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1 INTRODUCTION 

We do not know what fraction of stars form in dense star clus- 
ters. Part of the problem is that there is no agreed definitions of 
what 'dense' means and what a star cluster is (Portegies Zwart,| 
McMillan & Gieles 2010). In an attempt to shed light on this situa- 
tion Bresse rt et al.| ( [2010] > studied a sample of young stellar objects 
( YSOs) in the solar neighbourhood. They calculate the surface den- 
sity E around each YSO by finding the distance to the 7 th nearest 
nearest neighbour, d-j, such that E = 6/(7rd 2 ). They find that the 
distribution of E is roughly lognormal with a peak at about 22 pc -2 
and a dispersion of 0.85. Because YSOs are very young (of order 1 
Myr) they conclude that this distribution reflects the density distri- 
bution at the moment of star formation. They concluded that stars 
form in a broad and smooth spectrum of surface densities and that 
only a small fraction of the YSOs form in dense clusters. 

In this paper we argue that the observed peak in the surface 
density distribution of young stars (at around 20 stars per square 
parsec) is an expected outcome from a wide range of initial clus- 
tering configurations. In Section|2]we argue that for 'typical' clus- 
ter scales in star forming regions (iV ~ 100; Lada & Lada 2003) 
and young ages (about 1 Myr) of YSOs such a surface density rep- 
resents the outcome of dynamical evolution (i.e. two-body relax- 
ation) from a variety of plausible initial conditions. In Section [5] 
we flesh out this argument by considering the factors that broaden 
this distribution (the range of surface densities in a given cluster to- 
gether with a realistic spectrum of cluster membership number). In 
Section [4] we present a discussion and conclude that the observed 
surface density distribution is exactly what one expects if the ma- 



jority of stars are born in clusters; the situation is however highly 
degenerate so that it is not possible to use this distribution to place 
unique constraints on the initial conditions. 



2 THE PEAK OF THE CLUSTER SURFACE DENSITY 
DISTRIBUTION 

The evolution of self-gravitating clusters containing a realistic stel- 
lar mass spectrum is well understood. The massive stars segregate 
to the cluster core due to dynamical drag and initiate core collapse: 
in the case of clusters numbering a few hundred stars, this occurs 
after about 10 or 20 crossing times ( |Giersz & Heggie 1996). Sub- 
sequently the cluster expands due to the transfer of energy from 
binaries and stellar mass-loss in the core by two-body relaxation. 
This expansion approaches a self-similar state in which the clus- 
ter ex pands more or less hom ologously and with radius scaling as 
t 2 / 3 |Giersz & Heggie|l996^ . In this state the two-body relaxation 
timescale tends to a fixed multiple of the cluster age (HenonjT965}. 

Self-gravitating systems are scale free in the sense that, for a 
given distribution of stars, one can re-scale their spatial separations 
and not affect their qualitative evolution: in fact, scaling the spatial 
separations by a factor / simply involves a re-scaling of time by a 
factor / 3 ^ 2 . We illustrate this in Fig.[I]where we show the evolution 
of the half-mass radius for clusters of 256 stars and with a range 
of initial radii. It is evident that, although qualitatively identical, 
the initially compact cluster starts its self-similar expansion earlier. 
Notably, however, all three clusters are approaching the same self- 
similar track by an age of a few Myrs. 
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Figure 1. Evolution of the half-mass radii (r^) of clusters with N = 256 in 
direct N-body integrations with NBOD Y 6 ( Aarseth 2003 1. The stellar mass- 
function is afKroupa (2001) distribution between 0.1 Mq and IOOMq and 
the stars do not evolve internally. The initial positions and velocities were 
generated from a King (1966) model with Wq = 9. In total 200 simulations 
were done and the median result for j\ is shown for different initial rj, ■ The 
full lines shows the result for all the stars and the dashed line for only the 
bound stars. 



Another way to state this is that clusters numbering a hundred 
stars end up with a mean surface density of around < 100 pc -2 at 
an age of a few Myrs provided they start off with mean densities 
in excess of that. We therefore deduce that the peak of the surface 
density distribution observed by Bressert et al. is comparable to 
what one would expect from dynamical evolution of stellar clusters 
with ages and richness found in nearby star forming regions {Allen| 
|et al.||2007| |Gutermuth et al.||2009| >. We now turn to considering 
effects which account for the finite width of the distribution, such as 
the range of cluster membership and the range of surface densities 
within a given cluster. 



3 THE SHAPE OF THE OBSERVED DISTRIBUTION 

3.1 Surface density distribution in individual clusters 

We start by considering the distribution of surface densities in a 
single cluster. To facilitate a tractable analytical approach we use a 
|Plummer| ( |1911^ model to describe the surface (number) density of 
the cluster, which has a surface density profile of the form 



£(i?) = E 1 + 



I? 



(1) 



Here R is the distance to the centre, Ro is the scale radius and En 
is the central (and maximum) surface density. The latter depends 
on the other variables and the total number of stars in the cluster TV 
as En = N/(ttRo). The number of stars per logarithmic surface 
density unit is given by 



dn 
dlnE 



, dn 
J dR 



dR 



dT,{R) 



(2) 



Using dn/dR = 2n RY1 and equation |TJ we find 
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If we generalise equation l[TJ to any cored function with an asymp- 
totic surface density profile with E oc i?~ 7 then the power-law 
slope on the right-hand side of equation (3) becomes 1 — 7/2. 
For the Plummer model 7 = 4, hence we find a logarithmic 
slope of 1/2 in equation ([3}. A shallower profile of 7 = 3 would 
give a slope of 1/3. Note that below a surface density of about 
A -2 En the expected number of stars falls to < 1 and we thus have 
A -2 < E/Eo ^ 1. For a cluster with 100 stars there is thus a 
range of four orders of magnitude in densities. It is interesting to 
note at this point that the observed spread in densities in YSOs is 
(only) four order of magnitude (Bressert et al. 2010). This range is 
what one gets from a single cluster with A ~ 100 stars. In Fig. [2] 
we show the E distributions for stars in clusters with different A 
(dashed lines). The sum of the distributions of clusters with differ- 
ent A and En gives the overall E distribution and we will explain 
in section[3?2|how this is computed. 



3.2 Population of clusters 

To get the surface density distribution of all the stars in a population 
of clusters we need to know the number of clusters with A stars 
(i.e. the cluster initial mass function, CIMF). We adopt a power- 
law scaling for the probability that a cluster has between A and 
N + dN stars with N in the range Ai ow ^ N < A up 



4>ci(N) = AN~ 



(4) 



We can get the total E distribution by multiplying <fi c \(N) by the 
surface density distribution for the individual clusters (equation [2]l 
and integrating over all clusters, 
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To be able to do this, we need to relate the central surface density 
En to the number of stars in the cluster A. Following the argu- 
ments of Section 2 we use the size evolution of clusters of a given 
age in the self-similar regime of cluster expansion. This leads to 
a solution that is independent of what is assumed for the initial 
mass-radius relation. In this expansion phase the half-mass relax- 
ation time-scale, r r h, grows (roughly) linearly with time t (Henon 



1965 1. A constant ratio of r rh /t implies R oc N~ 1/3 t 2/3 jSpitzer 
& Hart|1971} . Since the cluster expands homologously, the central 
surface density is proportional to A/i? 2 and we can thus relate the 
central surface density En to the number of stars in the cluster and 
the age as 

Here £1 is a constant of proportionality that is discussed in sec- 
tion [33] These relations can now be used to perform a change of 
variables in equation {3}. The result of the integration over En is 

So, up 

(7) 
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where £n, up is the central surface density of the most massive 
cluster (i.e. En, up = Eii _4//3 Aup 3 , equation]^} and £n, m m is 
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the lower integration boundary. There are two regimes for which 
different forms for Eo, m in need to considered. For E smaller 
than the central surface density of the smallest N cluster (E < 
Eo,iow = Eit _4//3 A r 1 ^) all clusters contribute to the distribution 
(see Fig. [2}. In this regime E , m m = E ,i O w For E ,i ow < E < 
Eo, up we need to integrate only over the clusters that contribute to 
the overall contribution, that is, Eo.min = E. Using this we find for 
the surface density distribution of a population of clusters 



An 
dlnE 
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1 - E 
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low 1 

^0, low < S < 1, 



(8) 



where tildes denote surface densities normalised to Eo, up . 

The distribution has the behaviour of a single cluster (equa- 
tion^ at low densities (E ^ Eo.low)- In the regime above Eo,iow 
(but << So, up) the distribution is flat: this can be readily under- 
stood since the CIMF with —2 contributes equal numbers of stars 
in equal logarithmic bins of N (and hence, via equation 6, in equal 
logarithmic bins of E(i?)). The peak and width depend on the val- 
ues of Eo,iow and Eo, up , which both scale with Ei. So what re- 
mains to be done is to derive the constant of proportionately in Ei 
that sets these values. 



3.3 Constant of proportionately Ei 

To complete the model we need to know what the constant of 
proportionality Ei is in the relation between surface density and 
age (equation [6]l. Because we assume a functional form that goes 
through the origin we thereby also implicitly define the time-scale 
for clusters to reach the self-similar expansion phase. There are sev- 
eral factors that determine how fast clusters reach the asymptotic 
track, such as the density profile, the primordial binary star frac- 
tion, the amount of primordial mass segregation, the initial virial 
ratio, etc. Our aim in this section is the find a single value for the 
constant Ei that is consistent with different numerical results. 

We derive the constant of proportionality Ei from the results 
of two independent sets of iV-body simulations of clusters expand- 
ing under the influence of two-body relaxation. Note that these 
models neglect the presence of a gaseous component. This assump- 
tion is discussed in section [4] Let us start with the expression for 
the half-mass radius of an expanding cluster 



(9) 



/ , \ 2/3 

From the iV-body models shown in Fig. [T] we find that rh,i = 
1.3 pc. From direct iV-body integrations of star clusters with 
8 192 ^ N ^ 131 072, a full stellar mass-spectrum and includ- 
ing the effect of mass loss of the individual stars, it was found 
that expansion is important already at an age of a few Myrs for 
clusters with initial r r h less than 10 Myr ( [Gieles et al.||2010| >. 
From their Fig. 2 an expansion of rh/Vh(0) — 2 was found at 
1 Myr for a cluster with N = 8192 and an initial half-mass ra- 
dius of r n (0) ~ 0.086 pc (the initial r n was computed from the 
mass and initial half-mass density of the cluster). This implies 
r n ,i ~ 3.2 pc. The fact that this value is more than a factor of 2 
larger than the estimate from Fig.[T]illustrates that this parameter is 
quite uncertain. The difference is probably due to the fact that the 
N — 8192 has more massive stars compared to the N = 256 clus- 
ter, which speeds up the dynamical evolution (Gieles et al. 2010). 
|Moeckel et al.| p012) used direct iV-body integrations to evolve 
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Figure 2. Simple model for the distribution of surface densities around 
stars. All stars are in clusters that are expanding due to two-body relax- 
ation. The dashed line show the surface density of the individual clusters 
with different TV (equation^ weighted by the CIMF. The peak density of 
each cluster is given by equation {6} and the sum (full line) follows (equa- 
tions |8j from integrating over the cluster mass function, where we used a 
— 2 power-law distribution between 10 and 500 stars for the clusters (equa- 
tion]^. An age of 2 Myr was adopted. 



small TV systems (a few hundred stars) that were selected from the 
outcome of smoothed particle hydrodynamics (SPH) simulations 
of a star forming molecular cloud. From their Fig. 12 we derive 
1-3 % 7"h,i/pc < 4.5. The large variation in r n ,i might be because 
these clusters have undergone varying amounts of dynamical evolu- 
tion in the hydrodynamical part of the models. Also, these clusters 
start off with varying degrees of segregation of the massive stars 
towards the centre. 

Having compared the different results, we adopt r n ,i = 
2.5 pc. We are aware that there is an uncertainty of about a fac- 
tor of two in this value. For a Plummer model the scale radius Ro 
is about a factor 1.3 smaller than the half-mass radius. We can then 
write Ei = (7r.R0) -1 — 0.09 pc -2 . We use this and equation ^ 
to find Eo,iow and Eo jUp and to compute the total distribution of E 
using equation {8}. The result is shown in Fig. |2]for t = 2 Myr. It 
is not exactly log-normal, but instead asymmetric and skewed to- 
wards low densities. This is in fact something that was found in the 
observations of Bressert et al. (shown in Fig. [3}. 

The peak, or rather the flat top, of the distribution depends 
on the density of the average cluster in the CIMF (Fig. [2]( and the 
age of the cluster population. It shifts to lower densities as clusters 
expand. Therefore, the location of the peak in this model is insen- 
sitive to the details of the mass-size relation of clusters at birth or 
the densities around the stars when they form. This is only true if 
the evolving densities of all clusters have reached the asymptotic 
Eo(iV, i) relation of equation |6j, that is, all clusters must have 
been denser in the past. In this analytic model the peak is close to 
the 20 pc -2 peak of the observed distribution at an age of 2 Myr. 
In the next section we make a slightly improved comparison by 
generating the model in a Monte Carlo fashion. 
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3.4 Monte Carlo generated model 

Fig.[3]shows a direct comparison between a Monte Carlo generated 
sample of our model (points with error bars) to the result for the 
total of the three Spitzer surveys of young stellar objects (YSOs) 
(histogram) used in Bressert et al. (2010). All distributions are nor- 
malised to an area of unity. 

Each cluster population contains 10 clusters with the number 
of stars in the cluster randomly drawn from the mass-function of 
equation {4]l with N\ ow = 50 and 7V up = 500. This results in 
an average of roughly 10 3 stars in each cluster population, which 
is similar to the total number of sources in the individual surveys 
considered by Bresse rTet al.| J2010) . Note that the Orion Nebula 
Cluster (ONC) is excluded from their Orion sample and an upper 
limit of a few hundred stars is, therefore, probably appropriate for 
the Solar neighbourhood minus the ONC. For each star the sur- 
face density is determined by finding the distance to the 7^ nearest 
neighbour as was done for the observations (section[TJ. An age of 
2 Myr is adopted. We generate 1000 cluster populations and show 
the median values and the boundaries containing 67% of the points 
around the median as circles and error bars, respectively. The (blue) 
solid line shows the analytical result of section [3^2| 

The peak values of the simulated distributions are remarkably 
close to the observed peak. We also show the analytical model re- 
sults for ages of 1 Myr and 4 Myr, to show how the peak evolves 
in time. |Bressert et aL]{2010[ l find a small offset in the density dis- 
tributions of the class I and the class II objects, with the class II 
objects being slightly less dense. Because class II objects are older 
this could thus be interpreted as dynamical evolution (i.e. expan- 
sion) Q There are several physical effects we do not consider here, 
such as the effects of multiple overlapping clusters and the effects 
of the presence of gas on the cluster dynamics. Although it is not 
our aim to develop a model that reproduces the observations in de- 
tail we discuss what the effect of the presence of gas could be in 
the next section. 



4 DISCUSSION AND CONCLUSIONS 

We showed that the distribution of surface densities of YSOs pre- 
sented by Bressert et al. (2010), that was used to argue that only a 
small fraction of the stars form in dense clusters, can be reproduced 
by an extremely simple model in which all stars form in dense star 
clusters. In this model the location of the peak of the distribution is 
a measure of the age of the cluster population and the typical mem- 
bership number of young clusters. Our model requires the clusters 
to form with higher densities than they have at the present day, 
which is a constraint on the physics of the star formation process 
which can be tested with future facilities, such as ALMA. 

There are several physical effects that we have not included, 
such as the details of the spatial distribution of YSO on the sky and 
the presence of gas. Although it is beyond the scope of this work 
to present realistic models that include all these effects, it may be 
interesting to discuss the omission of gas in our modelling in a bit 
more detail. This is motivated by the observation that YSOs trace 
the molecular gas they form from and the total gas mass in star 
forming regions can be 10 or 20 times the total mass in YSOs (e.g. 
|Kennicutt & Evans|2012| >, with average gas surface (mass) densities 

1 Bressert et al. consider this but conclude, in fact, that the similarity be- 
tween the distribution of the class I and class II objects supports their as- 
sumption that the spatial distribution of YSOs is primordial. 
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Figure 3. Surface density distribution of YSOs (histogram) as found from 
the surveys used by Bressert et al. (20101. Results of the Monte Carlo sim- 
ulations using 50 ^ N ^ 500 are shown as open (blue) circles with error 
bars. The (blue) solid line shows the analytical model of section |T2] for an 
age of 2 Myr, where as the (red) dashed lines show the results for ages of 
1 and 4 Myr. To conform with Bressert et al.'s figure we plot dn/d log £ 
rather than log(dn/d log £) as in Figure[2] All distributions are normalised 
to unit area. 



three or four times the surface (mass) density of YSOs ( Gut ermuth| 
|et al.|2009| >. 

A star cluster embedded in an external potential (e.g. gas) will 
have a higher velocity dispersion, which slows down the two-body 
relaxation that causes the expansion we discuss here. This is be- 
cause r r h depends on the stellar velocity dispersion a and the stel- 



lar density p as r r h oc c /p l Spitzer & Hart||l971 i. To estimate 
the contribution to a of the gas we assume two spherically sym- 
metric distributions, with the same functional form for the density 
profile. Using rh and rh, g for the half-mass radii of the stellar and 
gaseous component, respectively, and M and M g for the corre- 
sponding total masses. The velocity dispersion of the stars can then 
be expressed in the stellar and gas parameters as ( Spit zer|1969) 



GM 



1 + 



M 



(10) 



The first term on the right-hand side is due to self-gravity of the 
stars and the second term is the contribution of the gas. If the gas 
and stellar distribution have the same half-mass radius then we find 
(X oc (M + M g ) / rh, which is what is usually assumed for models 
that consider the removal of natal gas from an embedded cluster. 
Introducing r\ — M s /M and p = rh, g /rh we can compare the 
relaxation time-scale of a gas embedded system to that of a system 
containing only stars (i.e. what is assumed here) 



r r h (stars and gas) = I 1 



/'- 



3/2 



T rh (stars). 



(11) 



From this we see that T r h of a gas embedded system is longer if a 
significant amount of gas is present (77 > 1) which has a compara- 
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ble half-mass radius (/i ~ 1). But the increase of r r h depends sen- 
sitively on the size of the gaseous system in which the stars evolve 
(a dependence), so the effect of an additional gas component 
becomes negligible if /j, > 3. Note that 77 / /x 3 is the ratio of the 
(volume) densities of the stars and the gas within their half-mass 
radii. This is not the same as the local (volume) density contrast. 
Estimates of 77 and (surface) density contrasts are available in some 
dense embedded clusters but it is not trivial to estimate the value 
of (j, from these observations. If fi ~ 1 the effect of dynamical 
expansion can be overestimated by an order of magnitude. On the 
other hand, if /j, > 3, our assumption of gas free cluster evolution 
is justified. 

From SPH simulations it was found that gas dominated star 
forming regions are in fact gas-poor on the scale of the (sink) par- 
ticles (i.e. fj, » 1, Moec kel et al.|20 l2 |. Argume nts as to why 
this might happen can be found in |Kruijssen et al.| ( [2012) . Simi- 
lar results were obtained from adaptive mesh refinement (AMR) 
simulations (Girichidis et al. 2012). The results of these numerical 
studies support our assumption of pure stellar dynamical evolution. 

We conclude that the distribution of surface densities can 
not be used as evidence that not all stars form in dense clusters. 
Similarly, the agreement between this model and the observations 
should not be construed as an argument that all stars necessarily 
form in clusters. Bres sert et al.| j2010) argue that a (roughly) log- 
normal density distribution is evidence against a scenario in which 
stars form in distinct 'clustered' and 'distributed' star formation 
modes. Their argument is that the distribution would be bi-modal 
(or multi-modal) if there were distinct modes. However, if we in- 
terpret the model in Fig.[3]as a 'clustered' star formation mode, we 
can add a 'distributed' mode with a density of several tens of YSOs 
per pc 2 and some dispersion and the total distribution would still 
be unimodal. Our results support the suggestion of Bressert et al. 
( |2010} that a (local) surface density threshold is not a useful tool to 
separate clusters from field stars, because in our model all stars are 
in clusters and there is a range of about four orders of magnitude in 
surface density. 
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